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With only a few exceptions, the numerical simulation of cosmic and laboratory hydromagnetic 
dynamos has been carried out in the framework of the differential equation method. However, the 
integral equation method is known to provide robust and accurate tools for the numerical solution 
of many problems in other fields of physics. The paper is intended to facilitate the use of integral 
equation solvers in dynamo theory. In concrete, the integral equation method is employed to solve 
the eigenvalue problem for a hydromagnetic dynamo model with an isotropic helical turbulence 
parameter a. For the case of spherical geometry, three examples of the function a(r) with steady 
and oscillatory solutions are considered. A convergence rate proportional to the inverse squared of 
the number of grid points is achieved. Based on this method, a convergence accelerating strategy is 
C*") ■ developed and the convergence rate is improved remarkably. Typically, quite accurate results can be 

obtained with a few tens of grid points. In order to demonstrate its suitability for the treatment of 
dynamos in other than spherical domains, the method is also applied to a 2 dynamos in rectangular 
boxes. The magnetic fields and the electric potentials for the first eigenvalues are visualized. 
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I. INTRODUCTION 



The hydromagnetic dynamo effect is the cause of the magnetic fields of planets, stars, and galaxies |21ll24| . In the 
last decades much progress has been made in the analytical and numerical understanding of magnetic field generation 
OA ' in cosmic bodies. Only recently, the homo gen eous dynamo effect has been validated experimentally in large liquid 
^ | sodium facilities in Riga and Karlsruhe (l3l Il4l If5l l25l l35j . 

The usual numerical method to simulate hydromagnetic dynamos is based on the differential equation method. 
In the case of kinematic dynamo models, for which the fluid velocity v is supposed to be given and unchanged, the 
I . relevant differential equation is the induction equation for the magnetic field B, 

m. ^ = V x(vxB) + AB, (1) 

4^ ' with hq and a denoting the permeability of the free space and the electrical conductivity of the fluid, respectively. 
Oh | Note that the magnetic field has to be source-free: 

g I V • B = . (2) 

Assuming that there are no external excitations of the magnetic field from outside the considered finite region, the 
boundary condition for the magnetic field reads 



B = 0(r" 3 ) as r^oo. (3) 

For a qualitative understanding of Eq. Q one should notice that the magnetic field evolution is governed by the 
competition between the diffusion and the advection of the field. Without the advection term the magnetic field would 
disappear within a typical decay time t& = Mo^/ 2 , with / being a typical length scale of the system. The advection can 
lead to an increase of B within a kinematic time ty. — I /v. If the kinematic time becomes smaller than the diffusion 
time, the net effect of the evolution can become positive, so that magnetic field self-excitation can start. Relating 
the diffusion time-scale to the kinematic time-scale we get a dimensionless number that governs the evolution of the 
magnetic field. This number is called the magnetic Reynolds number R m : 

Rm = Hocrlv . (4) 

Depending on the particular flow pattern, the values of the critical R m , at which self-excitation occurs, are in the 
range of 10 1 ...10 3 . 

For the more complicated case of dynamically consistent dynamo models one has to solve simultaneously the induc- 
tion equation for the magnetic field and the Navier-Stokes equation for the velocity, in which the back-reaction of the 
Lorentz forces on the flow has to be included. 

A considerable part of dynamo research has been devoted to magnetic field self-excitation in finite spherical bodies, 
such as the Sun or the Earth. Fortunately, for the spherical case the boundary conditions for the magnetic field can 
be formulated separately for every degree and order of the spherical harmonics, so that the treatment of the magnetic 
fields in the exterior can be avoided. 
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This pleasant situation changes drastically when dynamos in other than spherical domains are considered. Then 
the correct treatment of the boundary conditions becomes non-trivial. In particular, this problem arises in connection 
with galactic magnetic fields, and in simulations related to the recent dynamo experiments 0, 0, |2^, which 
are carried out in cylindrical vessels. There are three ways to circumvent this problem: 

• The correct non-local boundary conditions are replaced by simplified local boundary conditions, e.g. "vertical 
field conditions" or "pseudo- vacuum boundary conditions" [5|, l31| , demanding that the magnetic field has only 
a normal component at the boundary. This method is very cheap from the numerical point of view, but it is of 
course not correct. 

• The real dynamo body is virtually embedded into a larger sphere for which the well-known boundary conditions 
for every degree and order of the spherical harmonics can be used. The region between the real dynamo and 
the surface of the virtual sphere is thought to be filled by a medium with a lower conductivity than that of the 
dynamo domain. Scaling this artificial conductivity to lower and lower values one can look for the convergence 
of the results. This method was successfully employed for the simulation of the Karlsruhe dynamo experiment 
[28l l29l ] , where the dynamo module has an aspect ratio (ratio of height to radius) close to one. 

• The Laplace equation for the magnetic field is solved in the exterior of the dynamo domain, and the interior 
solution is matched to the exterior solution by using the correct boundary conditions. This method, which was 
used for the simulation of the Riga dynamo experiment |32l |. is correct but numerically expensive. 

This unsatisfactory situation concerning the handling of boundary conditions was our main motivation to reconsider 
the integral equation method to dynamos in finite domains [3j| . The formulation of this method for the case of steady 
dynamos, which is nothing other than the application of Biot-Savart's law to dynamos, can already be found in the 
book of Roberts Interestingly enough, in Roberts opinion (j^lj, P- 74) this formulation did "...not appear, in 

general, to be very useful." The integral equation method was used in a few previous papers 0, ITol ITU Il2j , in which 
the effect of boundaries was mostly discarded, however. The "velocity-current- formulation" by Meir and Schmidt |2^| 
was intended to circumvent the numerical treatment outside the region of interest. However, the numerical focus of 
this work laid more on coupled MHD problems with small magnetic Reynolds number than on dynamo problems. 

A concrete result of our recent paper |33| was the formulation of a system of one-dimensional integral equations for 
a dynamo model with a spherically symmetric, isotropic helical turbulence parameter a in a finite sphere, and the 
re-derivation of the solution found by Krause and Steenbeck [2!| for the special case of constant a. 

This system of integral equations for the case of spherically symmetric, isotropic a is also at the root of the first 
numerical examples considered in this paper. Our present goal is to study and optimize the performance of numerical 
schemes to solve the integral equations for dynamos of this sort. The restriction to spherically symmetric a has the 
advantage that the equations decouple for every degree and order of the spherical harmonics. That makes our method 
comparable to the corresponding integral equation method for the radial Schrodinger equation |iL ITil ITtI ] . From there, 
and from other applications of the integral equation method Q, H, 0, 0, |2tJ , it is well known that the linear systems 
arising from the discretization of integral equations are generally well-conditioned. We present the numerical results 
of an integral equation solver with a convergence rate proportional to the inverse squared of the number of grid points. 
We also show how the convergence can be improved drastically by using a convergence accelerating strategy. 

Whereas these examples for the case of spherical geometry illustrate the feasibility of the integral equation approach 
and its equivalence with the differential equation approach they do not demonstrate any particular improvement with 
respect to the latter. The main advantage of the integral equation approach, its suitability for the treatment of 
dynamos in arbitrary domains, is therefore exemplified by another example which would be very hard to deal within 
the differential equation approach. Again, we consider an a 2 dynamo, but restrict the electrically conducting and 
dynamo active domain to a rectangular box outside which we assume vacuum. For such " matchbox dynamos" , we 
compute the first eigenvalues and visualize the magnetic field and the electric potential structure. It is shown how 
the first three eigenvalues, which are different for the case of different side lengths of the box, converge for the case of 
a cubic box. 



II. BASICS 



In this section, we compile the necessary formulae which are at the root of our numerical investigation. For details 
of the derivation we refer to our previous paper |33[ . 

Basically, our considerations are restricted to the steady case, i.e., to dynamos with growth rate and frequency equal 
to zero. Quite generally plj . the electromotive force (emf) in turbulent flows of conducting fluids can be written in 
the form 



F = vxB + aB-/3VxB 



(5) 
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The first term in this equation is the usual emf induced in a fluid flowing with the mean velocity v under the influence 
of the magnetic field B. The second term, aB, represents the effect of a helical turbulence, with a characterizing the 
helical part of the turbulence that can be derived in the framework of mean-field magnetohydrodynamics [Ul • The 
concept of the a-effect plays a considerable role in dynamo theory. The term /3V x B reflects the decrease of the 
electrical conductivity due to turbulence. 

In the steady case, the current density j is given by 

j = (t(P - Vp) , (6) 

where a is the conductivity of the fluid and ip is the electrostatic potential. With these notations, the coupled system 
of integral equations for steady kinematic dynamos can be written in the following form [sH . l33| : 



^ , \ Mof f F(r') x (r — r') _ . una f ,,,,,, r — s' 

= ~t L Ir- rT dV - ~t J s X J^f dS 



(7) 



D S 

with [Iq being the magnetic permeability of the free space, n(s') denoting the outward directed unit vector at the 
boundary point s', and dS' denoting an area element at this point. D and S indicate integrations over the domain of 
the fluid and its surface, respectively. 

Note that in the case of infinite domains with constant conductivity the electric potential does not appear. In this 
case the integral equation system reduces to Eq. (J2J) without the boundary term. A wealth of numerical applications 
of this formulation for infinite domains can be found in the paper by Dobler and Radler . 

In the following sections we will illustrate the general approach (7J and |(HJ| with two different applications. 

III. SPHERICAL DYNAMOS 

In this section, a 2 dynamos in spherical domains will be considered. For this case a wealth of quasi-analytical and 
numerical results arc available from the differential equation approach which can be used for comparison. 

A. The system of radial integral equations 

As usual, we split the magnetic field into a poloidal and a toroidal part according to 

Bp = VxVx^rj , (9) 

B T = Vxf^rj . (10) 

In spherical geometry the defining scalars S and T and the electric potential are expanded in series of spherical 
harmonics Yi m : 



s(r t e,4>) = $>™M>WM) , (ii) 

l.rn 

T(r,e,4,) = 5> m (r)y im (M) , (12) 

l.rn 

v (r,e,<f>) = 5> im (r)y im (M) • (13) 

Lm 

For the special case that the only dynamo source is a spherically symmetric, isotropic a-effect it was shown |33l ] 
that the system of integral equations Q and (JHJ for the magnetic field and the electric potential can be transformed 
into the following system of integral equations for the expansion coefficients s; m (r) and ti m (r) of the defining scalars: 



•im(r) = ^ 



r' l+1 /V+ 1 



a{r')U m {r')dr' + j — a {r')t lm {r')dr' 



(14) 



4 



tim(r) = Hov 



J+l 



a(r)si m (r) - -j^j a(R) s im (R) 



I + 1 r l+1 



21 + 1 i? 2 '+! 



,i da(r') 



dr 



— si m (r ) dr 



l + l f r' 1 da(r') 



21 + 1 



dr 



— si m (r ) dr 



I 



J+l 



da(r') 



21 + 1 J r ' l+1 dr 



— si m (r )dr 



(15) 



This system of integral equations (|14|) and (| 1 51) is equivalent with the system of differential equations and boundary 
conditions 



hsir. 



1 
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1 



d 2 Sl m _ UJ_+ 1) 

dr 2 

d 2 t lm 1(1 + 1) 



-Sir 



dr 2 
1(1 + 1) 



tii. 



■ a(r)ti m , 

— ( a( r ) ^ Slm 
dr V dr 



a(r)sir, 



t lm (R) = R^^j[\ r=R +(l + l) Slm (R)=0 
dr 



(16) 

(17) 
(18) 



if we set the eigenvalue A; of the differential equation system equal to zero which corresponds to the steady case. 

Note that the effect of the electric potential at the boundary is already incorporated in Eq. I|15|) . The corresponding 
term ensures that the boundary conditions (|18fl of the differential equation system are automatically fulfilled in the 
integral equation method. 

B. Numerical implementation 

In this section, we first present a numerical method for the solution of the coupled integral equations l|14|) and i|15fl • 
Then, in order to improve the convergence and accuracy further, a convergence accelerating strategy is developed. 



1. The basic integral equation solver 
We start with the integral equation system l|14f) and p5|l . Setting x = r/R, and introducing the notations 

i 

G s (x,x') 



— j— j < x' < x 

x l 

J+l 



(19) 



X < x 1 < 1, 



i X 



G t (x,x') 



21 + 1 x l 



I 



J+i 



2l + l x ' l +^ 

Eqs. Ijl4(l and ljl5[) can be rewritten in the following form: 



< x' < x, 



x < x' < 1 



(20) 



Slm(x) 



fi aR 2 
21 + 1 



G s (x, x') a(x') ti m (x') dx' 



(21) 
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tlm(x) = HoV 



a{x)si m (x) - x l+1 a(l) sj m (l) 



i + i xl+1 [ x ,iM^ Slm{x , )dx , 



21 + 1 J dx 

o 

i 

+ f G t (x,x') da ^/ s lm {x')dx 



x ' 1 da ^ s , ( X >) dx > 
x l dx' Sl ^ x > dx 



(22) 



Substituting Eq. J2IJ into Eq. J23 yields 



where we use the definition 



wti m (x) = a{x) J G s (x,x)a(x)ti m (x)dx 
o 

i i 

+ ^2JT\ Xl+1 j j ^^r^ X '%rn(x'')G s (x', X '')dx''dx' 







1 1 



/ I (h(x,x')G s (x\x , ')a(x'% m (x'')^ d ^-dx''dx' 







-x + / G s (l,x )a(x )ti m (x )dx 

o 



x' da(x') 



y-G s {x' ,x")a{x")ti{x")dx"dx' . 







x l dx 



w = {2l + l)/{^a z C z R z ) 



(23) 



(24) 



with C denoting a scaling factor of the function a(x) according to the new definition a(x) — Ca(x). Therefore, the 
integral equation system l|21l) and l|22|l is reduced to the single integral equation (|23|) . 

For the numerical implementation of Eq. we decided to choose the classical extended trapezoidal rule. Of 

course, more sophisticated treatments of the integrals by Gaussian quadratures or Clenshaw-Curtis quadrature are as 
well possible. Despite its simplicity, the extended trapezoidal rule is chosen since it can be easily used as a starting 
point of a convergence accelerating strategy to be discussed in the next subsection. 

Choosing N equidistant grid points xi = iAx, with Ace = 1/N, and approximating the integrals by the extended 
trapezoidal rule according to 



N 

J f{x)dx « ^(/(^-i) + f&)) A * 



(25) 
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we obtain for the discretization of Eq. (|23[1 the following expression: 



Wtl m {Xi)= J2f=l {&(Xi)CjG s (Xi,Xj)ce(Xj) A x 

+ ^ CfcCj da ^ k ^ a{x j )Gt{xi,x k )G s {x k ,x j ){Ax) 2 
fc=i ' '' 

+ WJ~T x i ¥1 X) c fc c i4 ^J!* a(Xj)G a (x k , xj)(Ax) 2 

k=l 

-^^^^-aix^G^x^Axf 

k=l i X 

-0.5^^a(a; J ) Cj G s (x,,^)(Ax) 2 }t ;m ( a ; J ) , (26) 

where = l,i = 1, 2, • • ■ , AT — 1, cjv = 0.5. Eq. (|26|l may be written in the following matrix form: 

At = wt, (27) 



where 



with 



A — (a,ij)NxN , (28) 



fly = a(xi)cjG s (xi, Xj)a(xj) A x 

*a(x 
dx 



+ J~] CfcCj d "^ gfc ^ a(xj)G t (xi, x k )G s (x k) x j )(Ax)' 1 



fe=i 

+ i7^^ +1 ^ c fcCj4^^a(a; i )G s (a; fe ,a; i )(Ax) 2 
-i[ +1 a(1.0)c,G s (1.0,ij)a(i J )Ai 

_y-fk^*) 5( ) G ( )(Aa;) 2 

k—l 1 

-0.5 da ^ tiix^c-jCsix^x^iAx) 2 . (29) 

This eigenvalue problem can be solved numerically. First, the matrix A is reduced to the Hessenberg form, then 
the QR algorithm can be employed to obtain the eigenvalues w of the matrix A and hence those magnitudes C of 
the functions a(x) for which steady dynamos exist. It should be pointed out that, except for the particular case 
a = const, the matrix A is non-symmetric, hence the appearance of complex eigenvalues should be expected. 

In the following discussions, the method described in this subsection will be called the integral equation solver 
(IES). 

2. Convergence accelerating strategy 

The basic idea of the convergence accelerating strategy has been applied in the numerical solutions of various integral 
equations This strategy, which is actually based on the Romberg scheme for the numerical quadrature by a 

extended trapezoidal rule, appears in the literature under various notations as extrapolation method llflj or deferred 
approach to the limit 0- 

In this subsection, this convergence accelerating strategy will be adapted for our eigenvalue problem in order to 
improve the convergence and accuracy of the integral equation solver. 

It will be shown in subsection 4.2 that the convergence rate of the integral equation solver can reach ~ N~ 2 . This 
is also the theoretical error estimation of the calculated eigenvalues obtained from the extended trapezoidal rule with 
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the step size 1/N under the assumption that the kernel is sufficiently differentiable [?]. So if w is the exact eigenvalue 

of the integral equation l|23|) and wj, ' is the eigenvalue calculated by the integral equation solver, it is expected 
that 

w = w { ° ) +fiN~ 2 + 0(N- 4 ), (30) 
where /i is a constant. In Eq. (|3(J[1 . if doubling the grid number, we have 

w = vffi + \^N- 2 + 0(N' 4 ), (31) 

where Wg 1 ' denotes the eigenvalue calculated by using 2N grid points. From Eqs. (|30J) and J2IJ, we obtain 

w = wf ] +0(N- 4 ), (32) 

where — (4wq 1 ' ) — u>q°^)/3. Therefore, approximates the exact eigenvalue w with an error which is 0(iV~ 4 ). 

The above idea can be extended and a triangular array of entries wY is obtained. The array is generated from 
the first column of eigenvalues obtained from discretizing the integrals in the integral equation l|23[) by the extended 
trapezoidal rule with grid numbers 2^N, (j = 0, 1, 2, • • • ). 

j = o 4 0) 

, _(1) _(0) 
] = 1 w w\ 

. _(2) _(1) _(0) 
1=2 Wq w\ w 2 

Q _(3) _(2) _(1) _(0) 
J = O Wq Wl W 2 W 3 



The entry w£ is placed in the (j + l)th position of the (k + l)th column (j, k = 0, 1, 2, • • • ). In general, the entries 
in columns other than the first are obtained by the recurrence relation 

^)=(4^ + 1 1 )-«;« 1 )/(4 fc -l). (33) 

This idea discussed in this subsection will be examined by two examples in Section 4. 



C. Numerical examples 

In this section we will treat some example profiles a(r) by the developed integral equation solver. In order to 
validate the accuracy of the results they are compared with results known from other methods. 

Note that one has to distinguish between the eigenvalues A which appear in the differential equation system 
(jlSJI . and the values C as they result from the eigenvalues w of the steady integral equation system H26|l . The 
eigenvalues A of the differential equation system comprise as the real part the growth rate and as the imaginary part 
the frequency of the magnetic field mode. Both parts have a physical meaning. In contrast to that, the values C for 
the integral equation system give critical values for the intensity of a, which are only meaningful if they are real. A 
complex value for C has no physical meaning, it might only indicate the existence of a complex eigenvalue A in the 
vicinity of the real part of the critical value C . 

1. Known results 

The example profiles that will be considered are the following (see Fig. [I): 

1. a{x) = C , 

2. a{x) = Cx 2 , 

3. a{x) = C(-19.88 + 347.37a; 2 - 656. 7lx 3 + 335. 52x 4 ), 
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FIG. 1: The three considered examples of a(x). 
C(-19.88 + 347.37a; 2 - 656.71a: 3 + 335.52a; 4 ). 



Example 1: a(x) — C. Example 2: a(x) — Cx 2 . Example 3: a(x) = 
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FIG. 2: First example: a(x) = C. Growth rates Re(\) for the eigenvalues with I = l;n = 1...8. 



where C denotes the magnitude of the functions. 

The first example represents the well-known Krause-Steenbeck dynamo model, defined by a(x) — C . Its eigenvalues 
C are known to satisfy the relation J; +1 / 2 (C') = 0, with J;+i/2 denoting the Bessel functions of degree 1 + 1/2 |2lj . 
In the first row of Table [I] we give the first three eigenvalues C for I = 1 which we compute by the programme 
"Mathematica" . In order to validate later the results of our integral equation solver, the results are given with 14 
digits after the comma. 

Fig. [21 shows the first eight eigenvalues A for 1 = 1, depending on C, which are labeled by the radial wavenumber n. 
These curves result from a differential equation solver based on the shooting technique from Numerical Recipes |27| . 
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For a(x) = C the results of the differential equation solver code were shown to be equivalent with the exact value at 
least until 8 digits after the comma. For the remaining two examples we believe the accuracy of this code to be at 
least in the same range. For the Krause-Steenbeck dynamo model, the eigenvalues A are always real (Fig. [21 . 
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C 

FIG. 3: Second example: a(x) = Cx 2 . Growth rates for the eigenvalues with I = l;n = 1...8. The merging points of two 
neighbouring curves indicate a transition to a pair of complex conjugated eigenvalues (the imaginary part, the frequency, is not 
shown here). 

This situation changes for the second example function, a{x) — Cx 2 . Fig. [3] shows again the real parts of A for 
the first eight eigenvalues for 1 = 1. It is clearly visible that the curves of two neighbouring eigenvalues merge at 
certain points. At these points, two real eigenvalues A turn into a pair of complex conjugated eigenvalues (although 
the frequency is not shown in our plot). The corresponding critical values of C are shown in the first row of Table ITD 
Here we give only 8 digits after the comma which are well justified from the accuracy point of view. 

Example 3 has been chosen as it shows a very rich spectral structure, with complex eigenvalues A at the critical 
points were the growth rates are zero. Actually, this function provides the astrophysically interesting example of an 
oscillatory a 2 -dynamo, meaning that the eigenvalue with zero growth rate is oscillating whereas all the other growth 
rates are less than zero [34|. Fig. ^ shows the spectral structure, with its merging and splitting points where two real 
eigenvalues turn into a pair of complex conjugated eigenvalues, and vice versa. More details close to the critical point 
can be seen in figure [SJ The corresponding critical values of C are shown in the first raw of Table IIIII The value in 
parentheses gives the frequency (the imaginary part of A) at the critical value of C . Here we give only 4 digits after 
the comma, as we are less interested in the accuracy problem than in the problem of complex eigenvalues. 

2. Results of an integral equation solver 

Table shows, for the case of constant cv, the eigenvalues of C. The first row shows the eigenvalues C resulting from 
a solution of the equation J 3 / 2 (x) = by Mathematical the remaining rows show the results of the integral equation 
solver (IES) and of two variants of the accelerated strategy for different grid numbers N. The first variant, which we 
call AS1, corresponds to the choice k = 1 in Eq. (|33[1 . the second variant, AS2, corresponds to k = 2. The numbers 
in the first column are the number of grid points. 

Based on these values, Fig. [S]shows the relative error of the results for IES, AS1, and AS2. For the IES, the error 
decreases as ~ N~ 2 . The error is larger for the eigenvalues with higher n. This is no surprise as the eigenfunctions 
for higher n are more structured. For a grid number 512, we can obtain an relative accuracy of the order 10 -5 . 
The accuracy can be increased dramatically if the accelerating strategies AS1 and AS2 are employed. For AS1 the 
convergence is ~ N~ 4 , for AS2 it is ~ N~ 6 . For the latter, and a grid number 512, we obtain a remarkable accuracy 
between 10 -15 and 10~ 12 . Note that, in order to get accuracies better than 10 -9 , it was necessary to use the fourfold 
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FIG. 4: Third example: a(x) = C(-19.88 + 347.37a; 2 - 656.71a; 3 + 335.52a; 4 ). Growth rates for the eigenvalues with I = 1; n = 
1...8, and for I = 2; n = 1, I = 3; n = 1. At the merging points complex conjugated eigenvalues appear, at the splitting points 
two real eigenvalues re-appear. This is a real oscillatory dynamo because the eigenmode which becomes critical first has a 
complex eigenvalue. 




FIG. 5: Details of Fig. H The mode with I = 1, n = 1/2 crosses zero at C — where it is oscillatory. All other mode become 
critical for higher values of C. 



precision option of the FORTRAN compiler. 

Now let us consider the convergence for example 2, a(x) — Cx 2 (Table ITT1 and Fig. 0). The convergence rates for 
IES, AS1, and AS2 are again ~ A^~ 2 , ~ A^ 4 , and ~ A^ 6 , respectively. The errors are typically higher than for the 
case of constant a, which may have to do with the non-symmetry of the matrix to be inverted. For AS1 and AS2 we 
have skipped the last points in Figs. Qb and 0c as we do not have results from the DES with a higher precision. 

The results for example 3 are represented in Table IIIII Below the first row that shows the DES results for the 
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TABLE I: Eigenvalues for a(x) = C, I = 1. The first row gives the results of Mathematica. The second row shows the results 
obtained by the integral equation solver (IES) with consecutively doubling the grid number starting from 8. The third row 
represents the results obtained by the accelerating strategy one (AS1): (4 Wq — wTq )/3, j = 0, 1, ■ ■ • ,5. The fourth row 
represents the results obtained by the accelerating strategy two (AS2): (16 u7< i+1) -TtJ< j) )/15,j = 0,1,2,3,4. 
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.49323421983453 
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.72333448381441 
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.89531156477993 


16/32 
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.49339836694047 


7. 


.72512802101689 


10. 


.90353754660825 


32/64 
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.49340876255332 


7. 


.72524403674743 


10. 


.90408464071915 
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.72524759016372 
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.90408594539681 


16/32/64 
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7. 


.72525177112947 
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.90412111365987 
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4 


.49340945787289 


7. 
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10. 
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64/128/256 


4 


.49340945790849 


7. 


.72525183692173 


10. 


.90412165929768 


128/256/512 4 


.49340945790905 


7. 


.72525183693745 


10. 


.90412165942685 



critical value of C we give for two of the columns the value of the frequency that appears at this critical point. For 
those complex eigenvalues A the integral equation method cannot work properly because it is restricted to the steady, 
non-oscillatory case. However, it is interesting to observe in the rows below that the existence of complex eigenvalues 
of A is mirrored in the existence of complex eigenvalues of C . These complex values are unphysical; nevertheless their 
real part is not far from the correct real part, and the imaginary part indicates oscillatory behaviour. 

IV. MATCHBOX DYNAMOS 

In this section, we consider so-called "matchbox dynamos", i.e. dynamos in a rectangular box which is filled by the 
electrically conducting fluid and surrounded by vacuum. This problem allows us to illustrate how to discretize the 
system of the integral equations (7J and (JSJ) in the general case of non-spherical domains. 

A. Numerical implementation 

In this subsection, we develop a numerical scheme to solve Eqs. Q and iJSJ directly for the matchbox dynamo. In 
doing so, we have to cope with the singularities of the kernels in this integral equation system. Actually, many well 
developed and efficient analytical and numerical methods are available in the boundary element method 0, |26| to 
solve the singular integrals, even for integrals with strong singularities. For details of the treatment of the singularities 
in Eqs. Q and ©, we refer to the appendix. 

1. Numerical scheme 

In this subsection, the scheme developed in Section 3.1 is extended to solve the matchbox dynamo problem making 
use of the basic integral equations @ and JSJ- Our scheme is very similar to the so-called constant element method 
which is widely applied in the framework of the boundary element method j^- The difference is that in our scheme 
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FIG. 6: Relative error of the outcome of the integral equation solver for a(x) = C. The comparison is made with results from 
Mathematica. (a) Simple integral equation solver (IES). (b) Accelerating strategy 1 (AS1). (c) Accelerating strategy 2 (AS2). 



the trapezoidal rule is utilized for the discretization of the integrals over elements in order to improve the convergence 
and accuracy. 

Firstly, Eqs. Q and JSJ) are rewritten in the form 



~^(s) = i- J M,(s, v')B q {v')dV n q (s')G q (s, S 'Ms')dS', 



(34) 
(35) 



ID JS 

where the conventional Einstein's summation has been used. D denotes the matchbox, S is the surface of D, r 



13 



TABLE II: Eigenvalues for a(x) = Cx 2 . The first row gives the results of the differential equation solver (DES) . The next rows 
show the results obtained by the integral equation solver (IES) by consecutively doubling the grid number starting with 8. The 
third row group represents the results obtained by the average scheme one (AS1): (4 w ( j+1) - w ( j) )/3, {j = 0, 1, • • • , 5). The 
fourth row group represents the results obtained by the average scheme two (AS1): (16 w^ +V> — W^')/15, j = 0, 1, 2, 3, 4. 



DES 
IES 

8 

16 
32 
64 
128 
256 
512 
AS1 
8/16 
16/32 
32/64 
64/128 
128/256 
256/512 

AS2 
8/16/32 
16/32/64 
32/64/128 
64/128/256 
128/256/512 



11.46714098 17.15742615 30.20482435 35.92762083 49.02331308 



10.85648799 
11.30339950 
11.42577229 
11.45677452 
11.46454789 
11.46649262 
11.46697889 

11.45237001 
11.46656322 
11.46710860 
11.46713901 
11.46714086 
11.46714098 



16.67863125 
16.96226288 
17.10779976 
17.14499828 
17.15431823 
17.15664913 
17.15723190 

17.05680675 
17.15631206 
17.15739779 
17.15742488 
17.15742609 
17.15742616 



26.63454167 
27.81152422 
29.60723780 
30.05610549 
30.16768664 
30.19554253 
30.20250407 

28.20385174 
30.20580899 
30.20572806 
30.20488035 
30.20482783 
30.20482458 



42.49377786 
33.41913581 
35.27304970 
35.76666260 
35.88755487 
35.91761522 
35.92512011 

30.39425513 
35.89102100 
35.93120024 
35.92785229 
35.92763533 
35.92762175 



79.11841524 
41.42370422 
46.58751516 
48.42948706 
48.87567516 
48.98645261 
49.01410104 

28.85880054 
48.30878547 
49.04347770 
49.02440453 
49.02337843 
49.02331718 



11.46750944 17.16294574 30.33927281 36.25747206 49.60545113 

11.46714496 17.15747017 30.20572266 35.93387885 49.09245718 

11.46714104 17.15742668 30.20482384 35.92762909 49.02313298 

11.46714098 17.15742617 30.20482433 35.92762087 49.02331002 

11.46714098 17.15742617 30.20482436 35.92762084 49.02331310 



TABLE III: Eigenvalues for a(x) = C(- 19.88 + 347.37a; 2 - 656.71a; 3 + 335.52a; 4 ). The first row shows the critical value of 
C resulting from the differential equation solver (DES). The values in parentheses in the second row (only for n = 1|2 and 
n = 6|7) give the frequency at this point. The remaining rows give the outcomes of the integral equation solver (IES). The 
complex values have no precise physical meaning. However it is interesting that the real parts are not very far from the correct 
critical value of C. 

~ n=l— 2 n=3 n=4 n=5 n=6— 7 

TJES 1.0000 1.6788 1.9318 2.4544 3.2223 
(4.984) (13.538) 

IES 

10 0.9581—1.0695 1.2174 1.8313 3.3357 2.9734±0.2389 i 

30 1.1124±0.1011 i 1.6230 1.9264 2.2805 3.2575±0.0990 i 

100 1.1260±0.1073 i 1.6738 1.9313 2.4384 3.3019±0.1649 i 
300 1.1273±0.1079 i 1.6783 1.9317 2.4526 3.3214±0.1777 i 
1000 1.1274±0.1079 i 1.6788 1.9318 2.4542 3.3237±0.1790 i 
3000 1.1274±0.1079 i 1.6788 1.9318 2.4544 3.3238±0.1792 i 



(x,y,z) T , r' = (x\y',z') T . We use the notation 

K p<q (r,r') = -« p (r')G,(r ) r')+u i (r')G i (r,r')i M t^Gi(r ) r')o(r') ) (36) 

L p (r,s') = e pql n q (s')Gi(r,s') 7 (37) 

M q {s,r') = e piq Ui{r')G p {s,r') + a(r')G,(s, r'), (38) 

for p,q,i = 1, 2, 3, with the definition 

= G 2 (r,r') = G 3 (r,r') = J^. (39) 

As usual, 5 pq denotes the Kronecker symbol, and e pqi is the Levi-Civita symbol. 

The matchbox can be expressed as [0,a] x [0,6] x [0, c], where a, b and c are the lengths of the three sides of the 
matchbox. Now we divide this matchbox into (N — l) 3 equally sized small boxes Dijk(i,j, k = 1, 2, • • • , N — 1), which 
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FIG. 7: Relative error of the outcome of the integral equation solver for a(x) = Cx . The comparison is made with results 
of a differential equation solver, (a) Simple integral equation solver (IES). (b) Accelerating strategy 1 (AS1). (c) Accelerating 
strategy 2 (AS2). 



can be described as [xj,Xj + i] x [yj,yj + i] x [zk, Zk+i]- The lengths of the intervals [xj,Xj + i], [yj,yj + i] and [z/c,Zfc+i] 
are denoted as Ax, Ay and Az, respectively. The six faces of the matchbox are also discretized in a similar manner. 
For the faces z = and z = c, they are divided into the small rectangles [x,, Xj+i] x [yj, t/j+i] (i, j = 1, • • • ,N—1); for 
the faces y = and y = b, they are divided into [xj,Xj + i] x \zk, Zk+i] (i, k = 1, 2, • • • , N — 1); for the faces x = and 
x = a, they are divided into [yj, yj+i] x [z*,, Zfc+i] (j, k = 1, 2, • • • , AT — 1). In the following, wc denote a representative 
of these small rectangles as Slj(i s = 1,2, ■■■ ,6) which can be expressed as [x\ s i , x\ e i+1 ] x [x2j,x%j +1 ] • The lengths 
of [a^%, £i s i + i] and [x^-, a^'j+i] are represented as Axi and Ax2, respectively. 

The magnetic fields "sit" on the A^ 3 grid points of the volume (including the surface), whereas the potential "sits" 
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on the 6N 2 — 12N + 8 grid points of the surface. However, if the grid point is on an edge or a corner, we take it as 
two or three different grid points by considering it as located on different faces. Hence, wc have a total of 6N 2 electric 
potential degrees of freedom. 

With these definitions, Eqs. (|34[) and l|35l) become 



B P W = E / K M {v,v')B q {v')dV' 

^ i'j'k' JD i'j'k' 

i.i=l i'j' JS i'j' 

\vW = i E / M q ( S ,r')B q (r')dV 

i'j'k' jD i'j'k> 

~T E E f ix n q {s')G q {s, S ')y{s')dS'. 
77 i„i=i i'j' Js \f} 



(40) 



(41) 



For the integrals over Di'j>k', the application of the trapezoidal rule leads to 

E / K p , q (r,r')B q (r')dV * 

i>j>k> jD t'i'« 

^ c l 'C r c k 'K Ptq (r,r i , j , k/ )B q (r i , j , k/ )AxAyAz, (42) 

i'j'fc' 

E / M q ( S ,r')B q (r')dV' * 

i'j'k' J D i'j'k' 

^ Ci'Cj'Ck'MqisjTi'fk^B^Ti'j/k^AxAyAz, (43) 

where Ci' is defined as ci = 0.5, cat = 0.5, = 1.0, i' = 2,3, • • • ,N — 1, rvj//./ = (x;/ ,yj/ ,Xk') T ■ Similarly, for the 
integrals over S\fh, we have 



E E f^L^s'Ms^dS' 



6 

6 



E E / n,( B ')G,(s, B , Ms / )dS' 

i»i=l i'j' 



E E^ c i' n ?( s i ? j') G 3( s ' s ^-')^(4i') A ^ A ^. (45) 

i s i = l i'j' 



Substituting Eqs. I|42I45|I into Eqs. (f^TTf) and lj4"T|) and letting r = ly^s = s-°-, we obtain 

47T 



B p (r ljk ) = ^ d'Cj'Ck'Kp^rijk^i'j'k^Bqir^j'k^AxAyAz 



i'j'k' 
6 



E E^ / ^' ,i p( r «fc' s i'i')v 3 (s!/j/)Aa;iAa:2, (46) 



i.i=l i'j' 

4tt 



i'j'fc' 
1 6 

TZ E E c ^'^^'j') G ^ s ^ s ^')^^'j') Aa; i Aa; 2: (47) 
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where i, j, k, i', j' , k' = 1, 2, • • • ,N, and i s , i s ± = 1, 2, • • • , 6. 

Note that when r^j, belongs to Z}j<i'fe' , a weak singularity of the first integral of the right hand side of Eq. (|4(J|) 
occurs. We can employ the strategy discussed in the appendix to deal with such a singularity. For example, we can 
eliminate a small box [a^+i — gAa;, Xi+i] x [j/j+i — gAy,y J+ i] x [zfc+i — |Az,Z{ ;+ i] from D^fc when dealing with the 
weak singularity caused by setting r to r^/. in the integral 



/ K p , q {v,v')B q {v')dV. 



The overall effect of doing so is equivalent to setting K p , q (rijt.,Yi'j'k') to zero in Eq. when = ly^-//./. 

A similar technique can be applied to handle the singularity appearing in the second integral of the right hand side 
of Eq. 141|) . For example, we consider the following integral 



f n q ( S ')G q ( S ^,s'M S ')dS'. 



Since the point belongs to , it results in a singularity of this integral. We can define a small piece of surface as 
^Itj = [ x i,i+i ~ jAjci, x [x2.j+i — \Ax2,X2,j+i]- When proceeding the discretization, we just replace by 
Slj — S l / i: j and neglect the small piece S^. This is also equivalent to setting G 9 (s-°, s-fj,) to zero when s*j = sjfjv. 

Note that similar procedures can be employed to avoid the singularities of the other integrals in Eqs. (Q) and JSJ. 

Eqs. (|46[) and l|47l) can be rewritten in the matrix form 

X B = aWEX b - DX„), (48) 
0.5 X v = HX B - AX„ (49) 



where 



X B = (Bi(r m ),B 2 (r m ),- • • , J B 2 (r WA r A r), J B 3 (r ArA r A r)) T , (50) 
X v = ^(8i 1 ) J v(si2) s ---,V(s^)) T . (5 1 ) 

A((i a - 1)N 2 + 1)N + j, (i sl - 1)N 2 + (i' - 1)N + j') 
= ^ * cj, n q (s\fy) G q {B*,B\fy) Ax x Ax 2 , (52) 

H((i s - l)N 2 + 1)N + j, 3N 2 (i' - 1) + 3N(j' - 1) + 3(fc' - 1) + q) 

— a> Cf c fe / M q (s^,r ifj/k/ ) Ax Ay Az, (53) 

E(3N 2 (i - 1) + 3iV(j - 1) + 3(fc -l)+p, 
ZN 2 (i' - 1) + 37V (j' - 1) + 3(fc' - 1) + q) 

= c 4 < Cj/ c fc / If p , g (r iii: ,r<, .,-,*.,) Ax Ay Az, (54) 



D(3N 2 (i - 1) + 3AT(j - 1) + 3(fc - 1) + p, - l)iV 2 + (»' - l)iV + j') 
= — Cj/ Cj/ L p (ryfc,s^) A»i Aara, (55) 

with = 1,2,3, i,j,k,i',j',k' = 1, 2, • • • ,N and t, = 1,2, • • • , 6. 
From Eq. Q49[l. we obtain 

X„ = (0.51 + A)- 1 HX B . (56) 

For the inversion of the matrix 0.5 I + A some particular care is needed. Physically, the electric potential is defined 
only up to an additive constant, which implies that the matrix 0.5 1 + A is singular. This difficulty can be removed by 
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applying the deflation method 0, which is widely used, e.g., in the context of electro- and magnetoencephalography 
[20|. Actually, the matrix A can be replaced by 

Al=A+ 6A^ 11 ' (57) 

where we denote by Ii a quadratic matrix of the order 6N 2 x 6N 2 whose entries are all equal to one. Thus, Eq. l(5?)| 
becomes 

X ¥ , = (0.5I + Ai)- 1 HX fl . (58) 

Substituting Eq. (JSHJ) into Eq. gSJ yields 

-J-X B = (E-D(0.5I + Ai)- 1 H)X B , (59) 

where R m is the magnetic Reynolds number. This eigenvalue problem can be solved by the QR method, which gives 
the critical magnetic Reynolds numbers and the corresponding modes of the magnetic field and the electric potential. 
Although this numerical scheme is presented for the dynamo action in the matchbox, it can be easily extended to 
solve steady dynamo problems in other domains. 



B. Numerical examples 

In the following we will treat a 2 dynamos with a constant value a = C within a rectangular box. The most 
interesting situation is with vacuum in the exteriour of the box. Only for the cubic box we consider also the case with 
the exteriour space having the same conductivity as the interiour, in order to compare the results with the analytically 
known ones for spheres of comparable sizes. 



TABLE IV: The first eigenvalue for a dynamo with a(r) = C within a cubic box of sidelength 2 for the cases of conducting and 
insulating exteriour. The first two rows give, for the sake of comparison, the analytically known critical values for an enclosed 
sphere with radius 1, and for an enclosing sphere with radius y3. The remaining rows show the numerical results of the integral 
equation approach for different grid point numbers N in one direction. Note that there is a threefold degeneracy of the first 
eigenvalue due to the symmetry of the problem. 





Conducting 


Vacuum 


C sphere 


3.506 


4.493 


C sphere. 1 


2.431 


3.116 


IES 






5 


3.524 


4.254 


6 


3.292 


3.996 


7 


3.170 


3.866 


8 


3.098 


3.793 


9 


3.052 


3.750 


10 


3.021 


3.723 


12 


2.982 


3.694 


15 


2.952 


3.678 



In Table ITvl we show for the cubic box the first eigenvalues C in dependence on the grid number N in one direction 
(the total grid number is then A 3 ). To the best of our knowledge, there are no values available in the literature to 
compare our results with. However, there is at least a plausibility check for our results. Imagine two spheres, the first 
one, with radius 1, being embedded neatly into our cubic box, the second one, with radius enclosing the box. It 
should be expected that the eigenvalues for the cubic box are between those for the two spheres. As can be seen from 
Table ITvl this is indeed the case, both for the case of a conducting outer space and for an insulating exteriour. 

The convergence of the eigenvalue for increasing N is illustrated in Fig. |S| We have made a fit of the eigenvalue 
data to the free parameters /, g, and h in the function C(N) = f + gN 3h . The parameter h gives the convergence 
rate for increasing N, whereas the parameter / gives a reasonable estimate of the true eigenvalue. For the case of 
conducting exteriour this value is 2.919 as compared with 3.506 for the enclosed sphere and 2.431 for the enclosing 
sphere. For the case of vacuum the value is 3.656 as compared with 4.493 for the enclosed sphere and 3.116 for the 
enclosing sphere. Hence, the results are plausible in both cases. 
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FIG. 8: Convergence of the three-fold degenerated lowest eigenvalue for the cubic box with increasing grid points TV for the cases 
of conducting and vacuum exteriour. The fit curves are C(N) = 2.919 + 39. 5447V -2 ' 599 for the case of conducting exteriour and 
C{N) = 3.656 + 91.3017V -3 124 for the insulating exteriour case, respectively. N is the number of grid points in one direction, 
hence the total number of grid points is TV 3 . 



TABLE V: The three first eigenvalues for a 2 dynamos in a rectangular box, in dependence on the ratios of side lenghts. 

N Side lengths 1. EV 2. EV 3. EV 



8 


1.0:1.0:1.0 


3.793 


3.793 


3.793 


8 


1.0:1.0:0.8 


4.072 


4.128 


4.128 


8 


1.0:1.0:0.6 


4.524 


4.674 


4.674 


8 


1.0:1.0:0.4 


5.322 


5.515 


5.515 


8 


1.0:0.8:0.6 


4.878 


4.956 


4.500 


11 


1.0:0.8:0.6 


4.728 


4.898 


4.934 


8 


1.0:0.8:0.5 


5.235 


5.350 


5.436 



As for the convergence rate, the value -1.041 indicates a faster convergence than \ogN/N which was found by 
Dobler and Radler [| (we use N for their total number of grid points in order to distinguish it from our number N of 
grid points in one direction). This better convergence rate should be attributed to the use of the trapezoidal rule for 
the integration instead of the constant element method, what we have also confirmed by comparative computations 
with the latter method. 

If we replace the cube by rectangular boxes with different ratios of the side lengths we can see in Table El that the 
three-fold degeneration of the eigenvalues is lifted. 

For the particular case of a rectangular box with the sidelengths ratio 1.0:0.8:0.6 and a grid number N = 11 we have 
visualized the magnetic fields and electric potentials belonging to the three lowest eigenvalues (Fig. [[J}. The structure 
of the magnetic field with its typical mixture of poloidal and toroidal components is clearly visible. It becomes evident 
that the field with the dipole axis perpendicular to the largest box face has the lowest eigenvalue. 

It should be noted that we have also checked the divergence-free condition of the fields and the curl-free condition 
in the exteriour. Both conditions arc fulfilled by the integral equation method with a reasonable accuracy. 



V. CONCLUSIONS 



We have used the integral equation method to solve numerically the steady kinematic a 2 dynamo problem in finite 
domains. 

For spherical domains, our approach is similar to the integral equation method for the solution of the radial 
Schrodinger equation. In this case, with only some tens of grid points the method provides reasonable results for all 
three considered example profiles a(r). The error decreases with the inverse squared number of grid points. With 
the use of a convergence accelerating strategy, the accuracy and the convergence rate can be significantly improved. 
Interestingly, even oscillating solutions of the dynamo problem which cannot be reproduced by our steady method 
are at least mirrored by complex eigenvalues for the dynamo numbers whose real part is close to the correct critical 
value. 
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FIG. 9: Magnetic fields within the box (left) and electric potentials at the box boundary (right) belonging to the three lowest 
eigenvalues of the o? dynamo in a box with sidelengths ratio 1:0.8:0.6. The corresponding eigenvalues of C are 4.728 (top), 
4.898 (middle), and 4.934 (bottom). Evidently, the eigenmode with the dipole axis perpendicular to the largest box face (top) 
is most easily excitable. 
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The particular suitability of the method to handle dynamos in arbitrarily shaped domains was demonstrated by 
the treatment of an a 2 dynamo in rectangular boxes. 

In summary, the integral equation method seems to be an attractive tool for the treatment of hydromagnetic 
dynamo problems. The robustness and accuracy of the method encourages to generalize it to the unsteady case and 
to more complicated dynamo models. 
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VII. APPENDIX: THE HANDLING OF THE SINGULARITIES 



In this appendix, we present the techniques to handle the singularities in the integrals appearing in Eqs. Q and 
{J in general and for the matchbox in particular. For the integral 

F(r') x (r - r') 
d |r-r'|3 ' 

with r belonging to D, we split the domain D into two parts, one being a small subdomain D e which is usually defined 
as a ball of a radius e centered at r, the remaining being the region D — D e . Hence the integral can be decomposed 
according to 

F(r') x (r - r') ^, = f F(r') x (r - r') ^ 



'"' " I D-D '" '*'' " 



The first integral of the right hand side of this equation is a normal integral without any singularity. For the second 
integral, the introduction of the spherical coordinates system (p, 0, <fi) leads to 



f F(r') x (r — r') , e T 2lI 

/ j dV' = J Q e J™ J Q T F(p, 9, 4>) x (sin 9 cos (j>, sin 6 sin 0, cos ( 

J Dr. I 1 " _ r l 



— (It — | () ( |( | () X l//. I/. ' (Mil 1/ UIM.I. Mill/ MU «./.«. MM'' 1 

I J), 

sin9d<pd9dp (61) 



Assuming that the function F(p, 8, <f>) is finite we see that the right hand side of Eq. (|61|l vanishes in the limit e — > 0. 

Therefore, the considered singularity is a weak singularity. In order to avoid such a singularity in the numerical 
computation, we can just discretize the region D — D e instead of D. The error caused by this procedure can be made 
as small as desired by taking a small enough value of e. The same procedure can be applied to the first integral on 
the right hand side of Eq. JHJ in the case s = r'. 

For the integral 



/ ^(s>(s') . dS ', 
J s — s' 



appearing in Eq. JSJ), with s on the surface S, note that the unit vector n(s') tends to be perpendicular to the vector 
s — s' in the case that s — > s', that is, n(s') ■ (s — s') tends to be zero. If defining S e as a small surface of a size e 
including the point s, we obtain Q (p. 69): 

lim / ^(s')n(s') • S ~ S, 3 dS' -> 0. (62) 



,|3 

s — s' 1 



A similar strategy as mentioned above can be employed to avoid such a singularity by discretizing the surface S — S e 
instead of S. 
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Now we consider the integral 

r ^(s')n(s') x dS', (63) 

r — s' |c> 



r — s' 



with the point r sitting on the surface S. Denote S e as a small piece of the surface S, which satisfies |s' — r| < e. 
Then the integral (|63[1 can be written as 

J s ^( s ')n(s') x ^-^3 dS' = J s _ Sc ^(s')n(s') x ^pr dS' 

i^^Xs'jx^rfS'. (64) 

The first integral of the right hand side of this equation is a normal integral with no singularity. As for the second, 
defining another small disk of a small enough radius 77 centered at r, we have 

lim / ^(s')n(s') x , r ~ S „ dS' = lim lim / ¥>(s'Ws') x -^—^7 dS' 
|r-s'| 3 ^^°Jse-sJ |r-s'|3 

f r — s' 
= os(r)n(r) x lim lim / - — dS' , 

e^O r^0j Sr _ Sii |r-s'| 3 

where < i] < e. Introducing the local polar coordinates, dS' = pd9dp, leads to 

lim / , r ^ S ,, dS' = - lim / -dp [ (cos9,sm9) T d9 
"-^Js.-s* |r-s'| 3 1^0 J v p p 7 

,2tt 

= -limln-/ {cos 6, sin 6) T d6. (65) 
n^o 77 J 

Since the integrals cos 9 d9 and sin 9 d9 always vanish, we have 



Therefore, 



lim / 7— — ^dS' = 0. (66) 
^°Js t -s v \r-s'\ 3 

lim [ (p(s')n(s') x , r ~ 8 „, dS' = 0. (67) 
^°Js e |r-s'| 3 

This indicates that we can also discretize the surface S — S e instead of S in order to avoid such a singularity in the 
numerical computation. 

For more details of the handling of the various singularities, one may refer to j2(J (pp. 7-16). 
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